( 


' v 


COMPUTATION OF CANONICAL CORRELATION AND 
BEST PREDICTABLE ASPECT OF FUTURE FOR TIME SERIES 

Mohsen Pourahmadl 1 
University of California, Davis 2 

A.G. Miamee 3 
Hampton University 
Hampton, Virginia 


/ys~o S'/ 


ABSTRACT: The canonical correlation between the (infinite) past and future 
of a stationary time series is shown to be the limit of the canonical 
correlation between the (infinite) past and (finite) future, and 
computation of the latter is reduced to a (generalized) eigenvalue problem 
involving (finite) matrices. This provides a convenient and, essentially, 
finite- dimensional algorithm for computing canonical correlations and 
components of a time series. An upper bound is conjectured for the 
largest canonical correlation. 
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1 . Introduction 


For many practical and theoretical problems in time series analysis, 
cf. Akaike (1975), Tsay and Tiao (1985), Pourahmadi (1985), it is of 
interest to know or compute p the canonical or maximal correlation 
between the past 

P - Xt.i.Xt)', 

and the future 

F - [X t+1 ,X t+ 2,...]' 

of a stationary time series CX t } , and the corresponding canonical 
component or the best predictable aspect of future. Using the familiar 
ideas from multivariate analysis, this task, requires computation of 
eigenvalues, eigenvectors and inversion of infinite matrices or 

operators . 

For ARMA processes canonical correlations and components can be 
computed (exactly) by solving linear systems of algebraic equations, cf. 
Helson and Szego (1960) and Yaglom (1983). For general stationary 
processes, Jewell et. al. (1983) have given an algorithm for 

computing (approximating) the canonical correlations as the eigenvalues 
of an infinite- dimensional (Hankel) operator, in the spectral domain. 

In this paper, we provide a time domain algorithm for computing 
(approximating) canonical correlations of a (nondeterministic) 
stationary process, which requires only solving linear system(s) of 
algebraic equations, cf. Yaglom (1965). For instance, in our approach, 
computation of p m , the canonical correlation between the (infinite) past 
P and (finite) future 

F m ” t x t+l> • • • < x t+m3 ’ • l^nK 00 , 

requires solving one (generalized) eigenvalue problem for two mxm 
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matrices G m and r m , cf. Theorem 2.3. 

Our approach relies primarily on the Wold decomposition of a 
(nondeterministic) stationary process; this makes it possible to reduce 
the genuinely infinite-dimensional problem of computation of p m or p to 
an, essentially, finite-dimensional problem; in addition to its 
computational simplicity, this approach also provides a procedure for 
computing p even when it does not exist as an eigenvalue of an operator, 
cf. Jewell and Bloomfield (1983) and Jewell et. al. (1983). 

In section 2, we develop a procedure for computing the best 
predictable aspect of (finite) future; the main result is Theorem 2.3. 
This result along with a simple fact about geometry of Hilbert spaces are 
used in Section 3, to give an algorithm for computing p and the best 
predictable aspect of the entire future. This procedure is applied to 
the well-known models fitted to the sunspot numbers series; it turns out 
that even for m-4, /> m provides a good approximation for p. 

An interesting and yet open problem in this area is that of finding 

I f 2 

1- is an upper 

bound for p, where a' 2 is the interpolation error of a missing value 
based on the other values of the process. Throughout this paper, we have 
emphasized computation of the largest canonical correlation; other 
canonical correlations and components can be computed by following a 
standard procedure in multivariate analysis, cf. Theorem 2.3. 

2. Best Predictable Aspect of (Finite) Future 
For many practical and theoretical problems it is of interest to 
find the best predictable aspect of the future of a system; when the 
system is modelled by a stochastic process (X t ) , the problem of interest 
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can be restated as that of finding the best predictable linear functional 
of the future values of the form 

X - “ c r X t+r , (2.1) 

r-1 

where m<® and c^ c m are (necessarily) unknown; when m - <*>, (2.1) 

should be viewed as the limit in the mean of finite linear combination s . ^ 
In general, this is a hard problem to solve. 

For the time being, we deal with the simpler problem of finding the 
best linear predictor and prediction error of X in (2.1) when m and 
C]_,...,c m are known, and then in the next section we show how the 
solution of this apparently simpler problem can be employed to resolve 

the more difficult problem of finding the best predictable aspect of the 
future. The need for prediciton of linear functionals of the form (2.1), 

with known m and c^ c m , arises when the forecaster is interested not 

only in a forecast of Individual future values but also in forecast of a 
linear combination of m future values and a confidence interval for it. 
For example, if sales are recorded monthly, the forecaster might be 

interested in the forecast of next year's total sales (m - 12, c^ - ... - 

c 12 " 1) . or one might be interested in forecasting the average of some 

future values (c^ - ... - c m - 1/m), etc. 

Note that when c^ - ... - c m „i - 0 and c m - 1, then X - X c+m , and 
the prediction problem of X t+m can be solved in the time domain by using 
the Wold decomposition of [X t 3; in fact, with 

X t - 2 b k c t . k + V t , b Q - 1, Z b k 2 < ®, a 2 — 
z k -o K C K C ° k-0 k 


Var ( e t ) 


( 2 . 2 ) 
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representing the Wold decomposition of CX^}, where I s 

innovation process of CX t ) and CV t ) a deterministric process uncorrelated 
with Cc t ), the best linear predictor of X t+r is given by 


^t+r ” 2 I>k c t+r-k + V t+r , 

k-r 


(2.3) 


and its (mean square) prediction error is 

Var(x t+r - X t+r ) - a 2 ^ b’. (2.4) 

Note that (2.2) also gives rise to the following representation of 
r - [ 7 . .]. . , , the covariance matrix of CXt-}: 

T - a 2 TT* + T y (2.5) 

where T - (bj.jJ^ ^ with bj ** 0 for j < 0, and T y is the covariance 
matrix of the deterministic process CV t }. As it is expected, the 
prediction problem of the more general linear functional ( 2 . 1 ) also 
hinges on the Wold decomposition of CX t }. Indeed, from (2.1) and (2.2) 
we have 

® m » m 

^ l ^ 1 c r^ ) r+k-mJ c t+m-k + ^‘ 1 < rr^t+r* 

k-0 r-1 r-1 < - 



from this , X the best linear predictor of X based on P is 
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Consider the harmonizable process {X^, t c R) <=■ Lq(P) given by 
and its spectral bimeasure which is induced by Z, i.e,. 


e it6 Z(d0) 


F (A , B) - E Z(A) Z(B), 

2 

We claim that the corresponding spectral domain L (F) in this case is not complete. 
Verification. By our Eemma there exists a nonzero vector in H (-“) which does 

— y 

not have a series representation as in (A), Take one such vector V, Since V 

is clearly in H^(0) there exists a sequence L a£ Y_^ “V n of finite linear 

2 

combination of Y, 's; k 0 which converges to v in L n (P). We can write 

iC (J 

v - r f (0)z(de) 

n 'r 

where the nonzero functions f are defined on positive integers with f (k) = 

n n 

af. By our Theorem in section 2 we have 
k 


l f n " f m ' *F 


v - v 
n m 1 


Now since v converges to v and hence is Cauchy so is f , However this 
n n 

2 

particular sequence f^ of functions in L (F) does not converge to any element 
2 

f in L (F) . Because otherwise another application of the Theorem in section 


2 shows that f is in L 1 (Z) and 


H £ n - f ll F - 'I ■ l^n-j> 2 N • 


Thus yie see that V also converges to ] fdz , So 

n A. 


-f f dZ — ? f(i) Z ({!}) - Z f(i)\. Y ,, 
/R i-0 i-D 1 ~ 


which contradicts our choice ofV, 

REMARK 1 . Our example shows that the main result of [7] claiming the 
completeness of the spectral domain of any multivariate weakly harmonizable 
process is false even for a univariate strongly harmonizable process. 

REMARK 2 . We feel that the error in [7] occurs in lines 8 and 9 of the 
second column of page 4612, where the existence of a "certain projection onto 
a subspace" is asserted and a reference to page 33 of [90 is. made to 
support it. In view of the results established in this note the results in 
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which is a quadratic form whose matrix is the matrix of prediction 
errors. For computational purposes, it is important to note that the 
matrix G m is, indeed, the upper left mxm submatrix of the matrix 
G - </ 2 T'T, (2.8) 

where the (infinite) matrix T is as in (2.5), this provides a simple 
method of computing G m when the moving average parameters b] > ,b2,... are 
known or the task of Cholesky factorization of the covariance matrix T 
is accomplished. In the following r m also stands for the upper left mxm 
submatrix of T. 

The measure of (linear) predictability of any function X is usually 
defined as 


A(X) - 1 


Var(X-X) 
Var(X) ’ 


cf. Jewell and Bloomfield (1983). 

Next, we summarize some of the previous results. 

Lemma 2.1 . Let CX t ) be a nondeterministic stationary process with 

covariance function and moving average parameters b Q - 1, 

m 

bi,b 2 X - 2 c r x t+r w h®re m < ® and ci,...,c m are given real 

r-1 

A 

constants. Then, with X denoting the best linear predictor of X based 
on the infinite past X t ,X t . i , . . . , we have 

(a) Var(X-X) - c'G m c. 

(b) the measure of (linear) predictability of X is given by 

c'G m c 

A(x) " 1 * 7^1 • • 

Remark 2.2 . For m-1, the measure of predictability of X - X t+ ^ has the 
simple form A(X t+ i) - 1 a 2 /7o> cf- Lemma 2.1(b), since a 2 - 

expCjlog f (A) dA/27r3, where f(A) is the spectral density of the process, 


it follows that A(X t+ ]^) can be expressed in terms of the density of the 
process. However, for m>l , it seems difficult to find expressions for 
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A(X) in terms of the density. 

Next, we find the best predictable aspect of the future for a given 

m. In view of Lemma 2.1(b) this amounts to finding ci,...,c m such that 
m 

for X - 2 c r X t+r , A(X) is maximized. The next theorem shows how this 

r-1 

can be reduced to a standard (generalized) eigenvalue problem. 

Theorem 2.3 . Let CX t 3 be a nondeterministic stationary process and 
m 

X - 2 c r x t+r> for m>l fixed. Then X is the best predictable aspect 

r-1 

of the m future values X t+ i , . . . ,X t+m , if c - [c^ c m ] satisfies 

(G m - AT m )c - 0, for some A € R. (2.9) 

More precisely, let A^<. . .<A^, (k S m) be the distinct roots of the 
determinantal equation 

det(G m - AT m ) - 0 (2.10) 

and c(i) c (m) ^e c ^ e corresponding orthonormalized eigenvectors, i.e. 

c (i>r‘m c (j) “ ^i,j> ^ • j“l > 2 , . • • ,m. 

Then, c^i)F m . with F m - [X c+ ^ X c+m ] ' is t * ie best predictable 

aspect of future with the measure of predictability 

A(X (1) ) - 1 - A 1( 


and in general X^) - C(i)F m is the i c ^ best predictable aspect of future 
with measure of predictability given by 
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A(X (1) ) - 1 - Xi. 

Proof . Note that the problem of maximizing A(X) over the variation of c 
is equivalent to minimizing c'G m c subject to the side condition c'r m c - 1. 
Now, the results follow either from using the standard Lagrangian 
multiplier method, cf. Rao (1973, p. 583), or a method based on the 
Hilbert-Courant maximization Lemma, cf. Johnson and Wichern (1988, p. 
441). ■ 

For the purpose of computation it is important to note that roots 
of (2.10) are the same as the eigenvalues of the matrix S m G m SQ, where S m 
can be chosen to be either the inverse of the symmetric square root 
of F m or the inverse of the Cholesky factor of In the computation 

that follows we have used the latter . For a given time series data 
set Xi,...,Xx, the moving average parameters bi,b 2 ,... can be estimated 
either by fitting ARMA models to data or factorizing the estimated 
spectral density, cf. Jewell et al. (1983). 

2.4 Example . 

(a) For X t - + Cfl* Var(c t ) “ 1. m - 2, we have 

^2 “ [[ * ° 2 “ [[ 2 )’ d6t(G2 * A172) “ (1 ' A)(1 * 3A) 

with roots and corresponding vectors 


Al - 1/3, 

C(l) “ [* 2/76, 1/76]'. 

a 2 - 1 , 

c (2) - ( 0. V72] ', 


and best predictable aspects 

X (1) - -2/76 X t+1 + 1/76 X t+2 , A (X (1) ) - 2/3 
X (2) - 1/72 X t+2 , A(X (2) ) - 0. 

Note that X( 2 ) is actually uncorrelated with X t ,X t _i 
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(b) For X t - c t + .216 1 * *36 c t-2» Var(ct) " If we have 


i 

x l 


c (i) 


1 

.7291 

[.817, 

-.703, .363, 

-.258]' 

2 

.9419 

[-552, 

.71 , .061, 

• 21 l 

3 

1.000 

[ o , 

0 , 0 , 

.922] 

4 

1.000 

[ o , 

0 , .929, 

-.109]' 


(c) For X t 

- .216X c .i + 

.36X t .2 - e t > v ar(c 

t> - 1. 

m-4 , we have 

i 



c ( i ) 


i 

.810 

[ .818, 

-.573, 

0 , 0 ]' 

2 

.912 

[ .499, 

.736, 

0 , 0 ] 

3 

1.000 

l *. no , 

.409, 

-.511, .952] 

4 

1.000 

[ .343, 

-.096, 

.886, .306]' 


(d) For X t 


.36c t .i, Var( 

£ t > “ 1. 

m-4, we 

have 

i 

x t 


_ c (i) 



i 

.9756 

[ .976, 

.351, 

.126, 

•43 ]' 

2 

1.000 

[*.552, 

-.113, 

.975, 

.148] 

3 

1.000 

[-.009, 

-.001, 

-.019, 

.986]' 
.058] ' 

. 4 

1.000 

[-.214, 

.94 , 

.230, - 


It is interesting and important to note the pattern of zeros in the 


c^^'s in examples (b) and (c) , and explore their relationship with 
those in Tsay and Tiao (1985). 


3 . Canonical Correlation 

v 

It is well-known that there is a close connection between prediction 
(regression) problems and the concept of correlation; the root of this 
phenomenon can be traced to the following simple property of the geometry 
of Hilbert spaces: For N a subspace of a Hilbert space H, and X an 

X 

element of H with denoting its orthogonal projection onto N, we have 
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(X, P*) - (X, X-X + P*) - |x|* - (X, X-P*) 
- ||X||* - «X-pJ|a. 


(3.1) 


In this identity, the relationship between the correlation (angle) of X 
and P^, and their distance ||X-P N || 2 (prediction error) is rather self- 
evident; in the following we shall make deeper use of (an extension of) 
this identity in developing an algorithm for computing p m and p, the 
canonical correlation between P and F m , P and F respectively; which 
allows us to reduce a genuinely inf ini te- dimensional problem to a finite - 
dimensional problem; for this we need the following two useful lemmas 
which are not necessarily new and their proofs might be around in the 
literature. Due to the importance of these lemmas in our work, in 
Section 4 we provide proofs for these lemmas. 

Lemma 3.1 . Let M and N be any two subspaces of L 2 (fi) , the space of 
square integrable random variables. Then, 

sup| Corr(X,Y)| - 1 - inf |X-pg|> , 

S3 M-i 


furthermore with >?(M f N) denoting the above quantity, we have 
p(Mi,Ni) < p (M, N^) < p (M,N) , for any and Nj£N, that is p(-,*) is 

an increasing (set) function. 

Lemma 3.2 . Let [X t 3 be a stationary process, with P, F, F m , p and p m as 
before, and P m - [X c _ m+ i , . . . ,X t ] , m>l. Then, 


(a) p - p( P,F) - iim p m . 

m*^o 


(b) p *■ iim P (Pm* 


iim iim p(P m , 


Fn). 
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Next, we state and prove the main result of this section. It is 
instructive to compare the result in part (b) with Theorem 1 in Jewell 
and Bloomfield (1983). 

Theorem 3.3 . Let (X t ) be a nondeterministic stationary process with 

covariance function £ 7 ^} and moving average parameters b Q - 1 , b^,b 2 

Then, 

(a) p m , the first canonical correlation between the (infinite) past P 
and (finite) future Fm- is given by 

Pm ~ 7l-A lm - 

where A]^ m is the smallest root of the determinantal equation (2.10), 

(b) As m-*», p m t p, in fact, 

p - sup p m - 71 -inf Xi m , 
m m 

and the best predictable aspect of the (entire) future F is equal to 

iim , where X(^) is as in Theorem 2.3. 

m-**> 

Proof . Part (a) follows from Lemma 3.1 and Theorem 2.3, by taking N and 

M as the closed linear span of entries of P and respectively. 

(b) follows from Lemma 3.2(a) and Theorem 2.3. ■ 

Due to the importance of p in many situations, it is desirable to 

find accurate bounds for it, whenever it is not possible to compute its 

exact value. This problem has been studied by Jewell et. al. (1983) and 

some elementary upper bounds for p are given in terms of certain 

components of the spectral density of the process. A sharp lower bound for p 

x t+l 

can be obtained from Lemma 3.2(a) and 3.1 by taking X - — 1 e F; 

y-Yo 


p > J \ , cf. Remark 2.2. 


(3.2) 
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To show that this bound Is sharp, note that for an AR(1) process 

X t - aX t _i + c t , a 2 - 1, |a|<l, 

1 

we have 7 0 - ]_. a 2 and the bound J 1 - a 2 / 7 0 - |a| is attained by p. It 
is much harder to find a sharp upper bound for p; however, motivated by 
(3.2) we conjecture that when p<l , then 

p < 7 1 - a' 2 / 7 0 , (3.3) 

where a' 2 is the interpolation error of X t+ i based on CX S ; s + t+1). 
We note that the bound (3.3) is attained for the aforementioned 
AR( 1) process, since in this case, by using a result of Kolmogorov (1941), 
we have 

a' 2 - (J £' 1 (B)d0/2x)' 1 - 1, 

10 - 2 

where f(0) - | 1 - ae | is the spectral density of the AR(1) process, 
and 

J 1 * *' 2 /7o - 1*1 • 

A more solid motivation for the bound in (3.3) is the fact that p(-,-)» 
cf. Lemmas 3.1 and 3.2, is an increasing (set) function of its arguments; 
therefore, replacing N(P) by - sp(X s ; s * t+13, one arrives at a bound 
of the form 

y 1 - K a' 2 /7o . 

for p, where K is a constant. Thus, the conjecture amounts to showing 
that K - 1. 

Remark 3,4 . The canonical correlation between P and F(k) “ (^t+k> 

Xt+k+1 ■•••)> k^l fixed, denoted by p(k), can be also computed by the 
procedure developed in this paper. In fact, for any m > k, and taking c - 
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[0 0> c k+l c m ] one can P rove results similar to those in Sections 2 

and 3 for Pm-k( k ) • which is the largest correlation between P and 
[Xk+i , . . . ,X m ] . It is evident that /> m -k( k ) -* p(k) as m -*• ®, cf. Theorem 
3.3; in this case p m . k (k) “ J 1 - *l,m-k( k ) • where A l m . k (k) is actually 

the smallest root of 

det ( G m-k * * r m-k] ” 

where G^-k is the (m-k)x(m-k) matrix obtained from G m by deleting its 
first k rows and k columns. 

Example 3.5 . The well-known sunspot numbers series has been studied by 
many people and various models fitted to the data are given in Table l f 
cf. Jewell et al . (1983). We have calculated p 4 , that is the canonical 
correlation between P and F 4 , and the corresponding canonical component, 
using the method of Theorem 3.3, see Table 2. These results are very 
close to the results in Table 2 of Jewell et al. (1983) which contains the 
value p for these models; this suggests that the rate of convergence of p m 
to p must be rather fast. For model 2, pi± is far from p reported in 
Jewell et al. (1983), this difference persists even when m is large; it 
should be noted that model 2 represents a nonstationary process, and it 
might be that for such processes our approximation may not work well as 
far as compution of p is concerned. Despite this, the canonical component 
for model 2 is almost the same as that for model 1 . 
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Table 1 




Model 


Source 

1 

x t - 

1.34X t .! + 

. 65X t . 2 - e t 

Yule, Box- Jenkins 

2 

x t - 

1.62X t .x + 

X t -2 “ 

Yule 

3 

x t - 

1 • 3X C . ! + . 
.19X t .4 + 

54X c . 2 + .15X t . 3 
. 24X t . 5 - . 4X t . g - c t 

Bailey 

4 

Xt - 

1.57X t .i + 

1.02X c . 2 - .21X t _ 3 - c t 

Box-Jenkins 

5 

X t - 

1.42X t _i + 

•72X t . 2 - c t - -15^.]. 

Phadke and Wu 

6 

Xt - 

1 . 25X t _i + 

. 54X C . 2 * . 19X t . 3 - c t 

Morris , Schaerf 


Table 2 


Model 

2 

'4 

Canonical Component 

1 

.8566 

X t - .36X t+1 

2 

.99 

X t - .38X t+1 

3 

.8602 

X t - .268X t+1 - . 107X t+2 

4 

.9149 • 

X t - .474X t+1 - .082X t+2 

5 

.8476 

X t - .296X t+1 - ,044X t+2 

6 

.8676 

X t - . 409X C+1 + .126X t+2 


4. Proofs of the Lemmas 

In this section we provide proofs of Lemmas 3.1 and 3.2. 
Proof of Lemma 3.1. It is obvious that 


fCorr(X, P^) ; X 6 M} i CCorr(X.Y); X € M, Y 6 N3 , 
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and therefore. 


sug |Corr(X, P^) | £ sug |Corr(X,Y)|. 

Yen 


( 1 ) 


Also, for any X € M we have from (3.1) that 


| Corr(X, Y) | < |Corr(X, P^) | , for all Y e N, 


and thus , 


sug | Corr(X, Y) | < sug |Corr((X, Pj) | . 
YEN 


(2) 


Now, from (1) and (2), we get that 


sug |Corr(X,Y)| - sug |Corr(X, PjJ) | 
YEN 


(3) 


For any X E M we have 


Corr(X, P") 


(x.ffi 11^1 8 l p »« 

Ixll^l ' |x||r*| "W 


PI 


J i-h-pjl 


where Y 
have 

thus 


M 


€ M with || y|| — 1; furtehrmore , from (3.1) we 


Corr(Y , P*) 




{|corr(X, P*)|; X E m} - {J 1- «Y-pJJ|| 2 ; Y E M, ||Y|| - l} , 


or equivalently, 


sup |corr(X, P X ) | - fl - inf |x-pjl> 

X£M |5pi 


( 4 ) 


lb 


The desired result, now, follows from (3) and (4). I 

Proof of Lemma 3,2: The sequence Cp m ) is bounded and nondecreasing, 

thus it is convergent and, in fact, 

iim p m - sup p m . (1) 

m-*» m 

Also, since the linear span of F m is a subset of that of F, we have 
p m < p, for all m > 1, 

and therefore, 

iim p m < p. (2) 

m 

To establish equalitly in (2), note that for any two finite linear 

n m 

combinations X - Z a k x t-k> Y - Z buX t+ u, we have 
k-0 k-1 

| Corr(X, Y) | < p m < sup p m . (3) 

m 

By taking supremum of both sides of (3), over all X and Y as above, 
we arrive at 

p < sup p m . (4) 

m 


The desired result, now, follows from (1), (2) and (4). Proof of (b) is 
similar to (a) . 
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